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ABSTRACT 

Aims. To develop a fully Bayesian least squares deconvolution (LSD) that can he applied to the reliable detection of magnetic signals 
in noise-limited stellar spectropolarimetric observations using multiline techniques. 

Methods. We consider LSD under the Bayesian framework and we introduce a flexible Gaussian Process (GP) prior for the LSD 
profile. This prior allows the result to automatically adapt to the presence of signal. We exploit several linear algebra identities to 
accelerate the calculations. The final algorithm can deal with thousands of spectral lines in a few seconds. 

Results. We demonstrate the reliability of the method with synthetic experiments and we apply it to real spectropolarimetric ob¬ 
servations of magnetic stars. We are able to recover the magnetic signals using a small number of spectral lines, together with the 
uncertainty at each velocity bin. This allows the user to consider if the detected signal is reliable. The code to compute the Bayesian 
LSD profile is freely available. 
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1. Introduction 


The phenomena taking place on the atmospheres of solar-type 
stars is controlled by magnetic fields. One of the most straight¬ 
forward ways of detecting and measuring such magnetic fields 
is through spectropolarimetry. The observation of the intensity 
(Stokes I), linear polarization (Stokes Q and U) and circular po¬ 
larization (Stokes V) allows us to remotely infer the strength and 
orientation of magnetic fields. 

During the last decades, we have been able to detect mag¬ 
netic fields in cool stars. The first efforts ( |Donati et al.||1992] l 
detected the presence of polarimetric signals in the most ac¬ 
tive of them. In these cases, it was necessary to reach a very 
high polarimetric sensitivity to detect these signals in individual 
spectral lines. In general, tough, the polarization level in almost 
all cool stars is below our detection limits and we require post¬ 
processing techniques to enhance the detection rate. The first so¬ 
lution to this problem was the development of multiline tech¬ 
niques ( |Semel| 1989[ Semel & Li |1996| l, who averaged many 
photospheric spectral lines to improve the effective polarimet¬ 
ric sensitivity. This technique works under the assumption that 
noise is uncorrelated and adds up to zero when many spectral 
lines are summed. Meanwhile, the spectral line itself presents 
correlations that will make it appear above the noise after adding 
many spectral lines. This idea was later improved with the in¬ 
troduction of the least-squares deconvolution technique (LSD; 
|Donati et al.|1997| , which has facilitated the detection of polari¬ 
metric signals in a wide variety of stars and the generation of 
a set of interesting surveys for the detection of magnetic fields 
in different stellar types (e.g.,|Wade et al.|2012| |Marsden et al.| 
[20T4l [Martins et al.|2015| ). 

Recently, the LSD method has been critically analyzed by 
Kochukhov et al. pOlOjl, pointing out the difficulties on the in¬ 


terpretation of the final LSD profile in terms of a normal spectral 


line. In the recent years, several variations of the LSD method 
have been presented. These methods try to overcome some of 


the fundamental problems with LSD. Lor instance, Sennhauser 


et al. (20091 suggested to introduce an improvement to nonlin¬ 


early deal with the presence of blends in spectral lines, much 
close to what really happens during the formation of the spectral 
line. [Carroll et al.j P008| l and Martinez Gonzalez et al.] P008| l 
generalized LSD to account for the deviations from the assump¬ 
tion of a single common profile for all spectral lines. This was 
done by using the principal component analysis (PCA) decom¬ 
position of all the considered lines and truncating the decompo¬ 
sition to get rid of the noise. Semel et al. ( 2009| l introduced the 
idea of multiline Zeeman signatures as detectors of the presence 
of magnetic fields. Linally, pCochukhov et al.[ ( [2010| l introduced 
multiprofile LSD, a generalization of LSD which models each 
spectral line as the linear addition of several profiles, much in the 
line of the PCA approach. These authors also discuss a simple 
modification of LSD to regularize the solution using a first-order 
Tikhonov regularization, which favors smooth solutions. 

In this paper we treat the LSD from a Bayesian perspective 
and consider the use of a flexible prior for the expected LSD pro¬ 
file. This introduces a regularization that automatically adapts to 
the data. The method is fast and competitive in computing time 
with LSD. One of the main advantages is that, given that the 
method we propose is fully Bayesian, we are able to give sensi¬ 
ble error bars to the detection of LSD signatures. 


2. Bayesian LSD 

The LSD technique is based on assuming that the circular po¬ 
larization profiles of a set of selected spectral lines from a star 
correspond to a common basic Zeeman signature. The spectral 
shape of this signature does not change and only the amplitude 
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is modified from line to line. The validity of this approxima¬ 
tion is based on two assumptions. First, the line has to be in 
the so-called weak field regime, which happens when the mag¬ 
netic field is sufficiently weak so that the Zeeman splitting is 
small as compared with the Doppler width of the spectral line. 
Second, the line has to be weak, so that no saturation effects 
are detected. Under these simplifying assumptions, the radiative 
transfer equation can be solved analytically and it can be demon¬ 
strated that the Stokes V profile is proportional to the wavelength 
derivative of the intensity pr ofile (e.g.,|Landi DeglTnnocenti &| 
|Landolfi||20041 l. In this case, [Donati et al.| ( |1997| l demonstrated 
that it is possible to write the Stokes V profile for line i as: 


Vi(v) = a/Z(v), 


( 1 ) 


where v is the displacement with respect to the central wave¬ 
length of the line in velocity units and o',- is a proportionality 


constant that depends on the properties of the line. Donati et al. 
( |1997| l also showed that these constants are a, - Ajdigi, with A, 
the central wavelength of the line, dj is the depth of the line, and 
gi is the effective Lande factor of the line. 

With a model for all the lines at hand, we can write the gener¬ 
ative model. It states how a given observed spectral line is mod¬ 
eled in terms of some parameters. Let us assume that we observe 
Mines spectral lines, each one sampled at velocity positions. 
Then, we model our measurements with: 


Vi(vj) ^ aiZ(vj) + eij. ( 2 ) 

The previous expression indicates that our measured Stokes V 
profile at a given velocity position and for a certain spectral line 
is just a weight times the common signature at this very same 
velocity with some added noise, e/j. This noise can potentially 
be different for all velocity positions of all the spectral lines and 
even accommodate the more general case in which there is non¬ 
zero correlation between different pixels. Note that the genera¬ 
tive model can be written in matrix notation as: 

V = WZ H- e (3) 

where the matrix W is given by: 


ai 

0 

0 

0 ■■■ 

ai 

0 

0 

0 ■■■ 

03 

0 

0 

0 ■■■ 

0 

Ol 

0 

0 ■■■ 

0 

02 

0 

0 ■■■ 

0 

03 

0 

0 ■ ■■ 


In other words, the matrix elements are given by Wij — 

^jdjj mod J- 


2.1. Standard LSD 


It is very simple to obtain the standard LSD method in a 
Bayesian framework. To this end, let the observations V rep¬ 
resent the Stokes V profiles of line i at wavelength j ordered as 
V = [yi(vi),V’ 2 (vi),V’ 3 (''i),---,Vi(v’ 2 ),V 2 (v 2 ),...]^. Likewise, 
we define the Z vector as the sampled common signature, given 
by Z = [Z(vi),Z(v 2 ),.. .Y■ With these definitions, we can apply 
the Bayes theorem to obtain the probability distribution for the 
common signature given all the observed data: 


F(Z|V) = 


p(V|Z)p(Z) 

F(V) 
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(5) 


where p(Z|V) is the posterior distribution for the common sig¬ 
nature given the observations, /7(V|Z) is the likelihood or sam¬ 
pling distribution (that gives the probability of getting a set of ob¬ 
served spectral lines conditioned on a common signature), p(Z) 
is the prior distribution for the common signature (that encodes 
all the a-priori information we have for Z(v)). Finally, p(V) is the 
marginal likelihood or evidence that normalizes the posterior to 
make it a proper probability distribution. Given that the marginal 
likelihood does not depend on Z, it is usual to drop it. 

Even though the equations are valid for the general corre¬ 
lated case, we assume for simplicity that the noise is uncorre¬ 
lated and with variance cr]j for line i at velocity point j. Under 
this assumption, the likelihood can be written as: 

p(V|Z) = Af(V|WZ, diag(o-2)) (6) 

where Af(x|//, C) represents a multivariate Gaussian distribution 
for variable x with mean vector p and covariance matrix C. The 
vector cr^ contains all the variances ordered exactly like in the 
case of V. Note that if all pixels share the same noise variance, 
equal to cr\, the likelihood simplifies to 

p(V|Z) = M(V|WZ,cr2l). (7) 

If we assume that all Z(v) profiles are equally probable a- 
priori, the posterior probability for Z(v) is 

p(Z|V) cx 7V(V|WZ, diag(o-2)) (8) 

where the posterior distribution is now seen as a distribution with 
respect to Z. Eq. Q represents the full probabilistic LSD solu¬ 
tion. Given the special structure of the matrix W, it is very sim¬ 
ple to compute the maximum a-posteriori point estimate (the one 
that maximizes the posterior distribution) and one ends up with 
the known LSD result of|Donati et al.|(|1997ll (see also |Semel| 
let al.|2009l l: 


ZLSD(Vj) - 


lYcrJaiViivj) 


(9) 


2.2. Bayesian hierarchical model 

Erom the Bayesian point of view, the classical LSD method is 
just the maximum a-posteriori solution to the inference problem 
when flat priors are used for the common signature Z(v). In other 
words, the value of Z(v) at one velocity position does not depend 
at all on the value at different velocity positions. This assumption 
is not capturing some knowledge that we have from the physics 
of the problem. We know that there must be some correlation 
between different velocity positions because the Stokes V pro¬ 
file are a consequence of the combined solution of the radiative 
transfer equation and the statistical equilibrium equations on the 
stellar atmosphere. The emission and absorption properties of 
the atmosphere inherently introduce correlation between differ¬ 
ent velocity bins. One consequence of this is the typical smooth 
shape of the Z(v) profile, that we should somehow force in our 
solution. 

In the Bayesian framework, we impose this condition by in¬ 
troducing a suitable prior distribution. A very convenient way to 
do it is by using a Gaussian Process (GP; |Rasmussen & Williams] 
2005) l. Generally, a GP is a distribution over the subspace of 
functions that fulfill certain conditions. Prom the technical point 
of view, once a function is sampled at a finite number of points, a 
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Fig. 1. Examples of functions extracted from GP priors with a squared 
exponential (solid lines) and Matem with v = 3/2 covariance functions. 
In both cases, ctz = 1 and different values of Az are used. Note that Az 
acts as a scale over which the function varies. The functions extracted 
from the Matern covariance function have a rougher behavior than those 
of the squared exponential one. 


GP for the function Z(v) is defined by the following multivariate 
normal distribution: 


p{Z\az) = N(Z\ti,K(az)), 


( 10 ) 


where fi is the mean of the distribution, which is usually set 
to zero, and K{az) is the covariance matrix, that depends on a 
set of parameters az- GP are then defined by the specific co- 
variance matrix. Instead of writing the full matrix, it is custom¬ 
ary to define a function (the covariance function) that gives the 
value of the matrix elements. This covariance function has to 
fulfill certain conditions so that the ensuing covariance matrix is 
proper (for instance, the covariance matrix has to be symmetric). 
Widespread covariance functions are the squared exponential: 


Kijicrz,Az) = cr|exp 


(d 


•v/ 


24 


or the Matern class of covariance functions: 


,1-v 


Kij{crz,^z) - 

r(v) 




dz 




V^IVi ■ 


dz 


( 11 ) 


( 12 ) 


where Ky{x) is the modified Bessel function (see jAbramowitz &| 
|Stegun|1972| l and v is a real number controlling the roughness of 
the resulting functions. Matern covariance functions converge to 
the squared exponential covariance function when v —» 00 . Other 
covariance functions can be found in Rasmussen & Williams I 
( 2005| l. Fig. [T] shows samples from a squared exponential GP 
prior with (Jz - I and different values of Az- In the case of the 
squared exponential covariance function, one can see that the 
functions are smooth and the value of Az defines the rate of vari¬ 
ability. 

The two covariance functions described in the previous para¬ 


graph are stationary (Rasmussen & Williams 2005 1 ), because 
they only depend on the difference v, - Vj, so that they are trans¬ 
lation invariant. The experiments carried out in Sections and 
[^demonstrate that these covariance functions do a good work. 
Although we do not pursue this objective here, it is possible to 
work with non-stationary covariance functions that are more lo¬ 
calized and could possible explain better some features in some 
specific velocity bins (similar to the work of|Czekala et al.|2014|l. 


A GP prior introduces the parameters orz into the inference. 
These parameters of the prior are termed hyperparameters and 
have to be considered in the posterior distribution. This is the 
reason why the inference model is hierarchical. Therefore, the 
posterior distribution of the hierarchical model is now given by 


p(Z,az|V) 


piy\Z)p{Z\az)p{az) 

P(V) 


(13) 


where we have added a hyperprior, p(az), for the hyperparam¬ 
eters. Given that the hyperparameters are of no interest for ob¬ 
taining the common signature, they should be marginalized (in¬ 
tegrated) from the posterior, so that we obtain the posterior dis¬ 
tribution for Z: 


P(Z|V) = 


J dazp(VlZ)p(Zlaz)p(az). 


(14) 


This integral cannot be computed in closed form in general. For 
this reason, we pursue here a Type-II maximum likelihood (ML) 
solution |Bishop| ( [T996| l, also known as empirical Bayes or maxi¬ 
mum marginal likelihood. This method is used for the analytical 
treatment of hierarchical models in which the posterior distri¬ 
bution is intractable. The strategy is to compute the marginal 
posterior for the hyperparameters: 


p(t^zlV) - -^p(az) f dZp(Y\Z)p(Z\az), (15) 

P(V) J 

and assume that this distribution is strongly peaked. Under this 
condition, one can obtain the set of hyperparameters az that 
maximize the marginal posterior p(azlV) and use a Dirac delta 
function as prior distribution: 


p(az) = d(az - az)- (16) 

In such a case, the integral in Eq. ( [I4| ) can be trivially carried 
out, so that we find 

p(Z|V, az) = ^p(V|Z)p(Z|az), (17) 

F(V) 

where we make explicit in the posterior for Z that we are condi¬ 
tioning on a specific value of the hyperparameters. 

The advantage of following this Type-II ML approach is that 
the solution becomes analytical. Making use of the definition of 
the likelihood of Eq. ([^ and the prior of Eq. ( [T0| i with p - 0 
and applying standard properties of the multivariate Gaussian 
distribution, we find: 


p(Z\az,y) ^ N(Z\pz,^z) 

where 


£7 = 


K(q'z)"' -HW^diagl^jW 


1 


Atz = EzWMiag ^ V 




Likewise, the marginal posterior for az is: 

p(az|V)=/7(az)W(V|//„,i:,), 

with 


(18) 


(19) 


( 20 ) 


= diagio-^) + WK(trz)W^ (21) 
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Fig. 2. Examples of two different shapes of Z(v) and how they are recovered using LSD (green curve) and Bayesian LSD (mean showed in blue and 
credibility intervals showed as blue shaded regions). The original true Z(v) is shown in red. The experiments use the Matem covariance function 
with V = 3/2. 


The objective now is to obtain the value of az that maximize Eq. 

This is easier if one takes logarithms, ending up with the 
maximization of 

logp(azlV) = logp(az) - ^ log|EJ , (22) 


again of size Nv x N^. These matrices are of a much smaller size 
than Eq, and the inversion is much less time consuming. Like¬ 
wise, using the matrix determinant lemma, we find: 


where |Eq| is the determinant of the matrix. 

From a technical point of view, the calculation of the inverse 
and the determinant is computationally heavy. The reason is that 
matrix Eq is of size MinesA^v x which is very large in 

practical applications. However, it is possible to efficiently accel¬ 


erate the computation using the Woodbury matrix identity (Pe- 
[tersen & Pedersen|2012| l, so that 

= diag(^ 


diag 


W 


K(ofz)^' -HW^diag — W 


-1 


W^diag, , 
(23) 


In this case, it is only necessary to obtain the inverse of the 
Nv X Nv matrix K(Qrz) and of the matrix inside brackets, which is 


log |Eq| = log 


K(Qrz)-' 


-H W^diag 



W 


+ log |K(Qrz)| + log |diag(0-2)1. 


(24) 


The inverse of the matrix inside brackets can be efficiently com¬ 
puted using the Cholesky decomposition which also gives as the 
logarithm of the determinant as a subproduct. 

As a final remark, note that many of the terms in Eq. (j2^ 
do not change during the optimization with respect to the hy¬ 
perparameters. They can be precomputed at the beginning of the 
optimization and only the term inside brackets in Eq. ( |2^ needs 
to be recomputed when the hyperparameters are changed. 
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Fig. 3. Results of applying the LSD (green curve) and Bayesian LSD (blue curve showing the mean and blue regions showing the credibility 
intervals) to observations of II Peg, both for Stokes I (left panel) and V (right panel). The experiments use the Matem covariance function with 
V = 312. The green and blue curves overlap almost exactly in the left panel when increasing the number of considered spectral lines. 


3. Synthetic experiments 

The first test of the method is done with synthetic data. We gen¬ 
erate two artificial Z(v) profiles, which are given by: 

Zi(v) = exp[-v^] 

Z 2 (v) = exp [-v^] -H ^ exp [-(v - 2.5)^]. (25) 

For convenience, both profiles are normalized to 1 /2 amplitude, 
multiplied by randomly generated a coefficients (with zero mean 
and unit variance) and then noise with diagonal covariance with 
variance cr^ = 1 is added to all the synthetic observations. A 
maximum of 200 synthetic spectral lines are used in both the 
standard LSD and the Bayesian LSD algorithms. For a single 
line, the signal-to-noise ratio equals one S/N - 1 /2. It is ex¬ 
pected that, in the presence of noise, adding Mines spectral lines 
would increase the S/N in proportion to VMines- Figurej^shows 
how both estimations approach the true common signature (rep¬ 
resented in red). The Bayesian LSD result is using a Matem co- 
variance function with v = 3/2. The standard LSD method gives 
the green estimation. It can be seen that, given that no correlation 


is imposed between velocity bins, a large oscillation produced by 
noise remains even when adding 200 lines. 


The mean of the marginal posterior for the common signa¬ 
ture, given by in EO- is represented in blue. One of the 
obvious advantages of working in a Bayesian framework is that 
the solution is given in terms of a probability distribution. This 
allows us to compute the probability distribution for each veloc¬ 
ity point. They are obtaining by marginalizing out all variables in 
Eq. (fTSll except for the value of Z(vj) at the velocity bin j of in- 
teresLUsing the special properties of the Gaussian distribution, 
these distributions are just one-dimensional Gaussian distribu¬ 
tions with variance equal to the (j, j) element of matrix We 
represent the distribution at each velocity point using the 68%, 
95% and 99% credibility intervals with blue shaded regions. 


Several properties of the Bayesian LSD solution are worth 
mentioning. First, the mean of the solution is very smooth. This 
is a consequence of the strong regularization introduced by the 
GP prior with squared exponential covariance function. Other 
covariance functions that give a less smooth solution can also 
be used. Second, the covariance function is stationary (it only 
depends on v - v') so that we find some variability in the inferred 
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Fig. 4. Same as Fig.j^for HR1099. 


solution when approaching the continuum. The reason is that 
the scale of variation Az is the same for all the velocity bins. 
Many other covariance functions can be used to diminish this 
behavior. Third, the true solution is always inside the credibility 
intervals. Fourth, the marginal posterior of Eq. ( fT4| ) can be used 
as a suitable distribution to carry out the modeling of the LSD 
profile. Using this distribution instead of the usual least squares 
solution (which assumes a diagonal covariance) makes the fitting 
much more robust. 


4. Real data 

The previous section demonstrates that Bayesian LSD can re¬ 
cover a much more robust common signature from the observa¬ 
tions, even when it is quite complex. In the following, we ap¬ 
ply the method to real observations. We consider three different 
cases, that are representative of what one can encounter in typ¬ 
ical observations. The case of II Peg is representative of a very 
strong signal that can even be detected in individual lines. The 
case of HR1099 is representative of a very complex common 
signature. Linally, 18 Sco is representative of a very weak signal 
in which, even taking into account thousands of lines, the final 
signal is very close to the noise level. Public spectropolarimetric 
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data for these three targets were downloaded from the PolarBase 
data base ( [Petit et al.|2014] l. 


4.1. II Peg 


II Peg is a very active RS CVn binary system, where a system of 
starspots has been mapped on the KlIV primary component. The 
spot system usually covers up to 40% of the surface ( Berdyuginaj 
et al. 19981. The observations have been carried out in July 2005 
with the ESPaDOnS spectropolarimeter ( jPonati et al. [2006) 1. Us¬ 
ing a total integration time of one hour, the peak signal-to-noise 
ratio of the original Stokes V spectrum (achieved at wavelengths 
around 780 nm) is close to 1,000 per wavelength bin. In agree¬ 
ment with the vast majority of published studies based on ES¬ 
PaDOnS data, the S IN used here is the one computed from the 
Stokes I spectrum obtained by adding up the four polarimetric 
exposures combined to produce the Stokes V spectrum. Eor such 
extremely active target, this noise level is sufficiently low to de¬ 
tect Zeeman signatures in a number of individual spectral lines 
formed in the photosphere and chromosphere ( [Petit|200'6l l. Here, 
we use a list of about 7,400 spectral lines computed from an 
LTE atmospheric model matching the spectral type of the pri- 
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Fig. 5. Same as Fig.j^for 18 Sco. 


mary component of II Peg ( |Kurucz|I993 i, keeping only photo- 
spheric lines deeper than 40% of the continuum level. 


Fig. shows the improvement in the detection of the com¬ 
mon signature Z(v) when the number of lines that are considered 
increases. The left panel shows the common signature in Stokes 
I, while the right panel displays the results for Stokes V. The 
lines that we have considered have been ordered on decreasing 
value of the weights a,-. Therefore, the panels with only 7 lines 
are obtained considering the lines that produce the largest Zee- 
man signal. The LSD profile is shown in green, while the mean 
of the Bayesian LSD detection is shown in blue, with the cred¬ 
ibility intervals at 68%, 95% and 99% levels displayed as blue 
shades of different grades. For the case of II Peg, the signature 
is detected at a very high signal-to-noise ratio with only 7 lines 
(0.1% of the total number of lines considered in the line list). 
The Bayesian LSD approach leads to a much smoother profile 
than the LSD one. The presence of credibility intervals is of rel¬ 
evance for discarding features like the one found at v =» -85km 
s ', which is compatible with zero in the Bayesian LSD case. 
Increasing the number of spectral lines produce a convergence 
of the Bayesian LSD and the standard LSD profiles, which is 
specially relevant in Stokes V. 


4.2. HR 1099 


Similarly to II Peg, HR 1099 is one of the most active RS CVn 
systems, with both components of the binary association clearly 
visible in high resolution spectra (while the primary alone is seen 
in observations of II Peg). The projected rot ational velocit y has 
been measured to be v sin i = 41 +0.5 km s * ( Donati 1999 1 . The 
magnetic structure of its active subgiant has been investigated 
in detail by [Petit et aT] ( |2004| l using Zeeman Doppler Imaging, 
showing the presence of large and axisymmetric regions with the 
magnetic field being azimuthal. Additionally, [Petit et al. ( 2004[ l 
also found clues for the presence of a differential rotation, much 
weaker than in the case of the Sun. The observation considered 
here was taken with ESPaDONS in August 2005. 


The line list used here is the same as the one previously em¬ 
ployed for II Peg, providing us with a good match to the photo- 
spheric characteristics of the evolved primary component, while 
the secondary is a much less active G5 main sequence star. 


The presence of several magnetized regions that also change 
with time produce Zeeman signals of large complexity. Figure 
1^ shows the inferred common signature for HR 1099, both for 
Stokes I (left panel) and for Stokes V (right panel). Several mag¬ 
netic components are detected around -65 (which is roughly the 
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radial velocity of the secondary star at that time), 20 and 40 km 
s“' (originating from the surface of the primary star of the sys¬ 
tem). These detections in Stokes V are also consistent with what 
is found in Stokes I. When only the 7 most magnetically sensi¬ 
tive lines are considered, the mean of the Bayesian LSD profile 
for Stokes V is very smooth, with the dashed regions marking 
the large uncertainty. A large part of the noise that is intrinsic to 
the standard LSD profile is absent. With only 74 lines, the fea¬ 
ture at positive relative velocities clearly appears. The quality of 
the detection increases drastically when considering more than 
10% of the lines of the line list. However, it is interesting to note 
that the uncertainty in the continuum regions does not decrease 
much when using 10% or 100% of the lines (third and fourth 
panels). In the presence of only pure random errors with constant 
variance, this uncertainty should be reduced roughly in propor¬ 
tion to VMines- However, the assumption of constant variance in 
these observations is not correct because the S /N changes much 
throughout the ESPaDOnS spectrum, being much lower at the 
borders of individual spectral orders and also near the bound¬ 
aries of the spectral domain of the instrument. This leads to a 
slower decrease of the final uncertainty. It is important to point 
out that our Bayesian LSD method is able to take this into ac¬ 
count and the ensuing final uncertainties do not decrease in such 
proportion. The immediate consequence of this is that we could 
be missing very weak signals due to the presence of systematic 
errors. 


4.3. 18 SCO 


18 Sco (HD 146233) is known to be one of the best solar twins 
among bright stars ( Melendez & Ramlrez|2007) l. In spite of its 
relatively low activity level, its surface magnetic field was de¬ 
tected by [Petit et al.| ( [2008] ) using high S/N observations. The 
Stokes V spectrum analyzed here was obtained with the NAR- 
VAL spectropolarimeter ( |Auriere| 2003| l in July 2008. This ob¬ 
servation is extracted from the time-series presented by |Petit| 
|et al.| ( |2008| l. Our line list is also identical to the one used in 
the same study. The LTE model from |Kurucz| ( |1993| ) was com¬ 
puted for a G2 spectral type, leading to a total of about 3,700 
photospheric spectral lines deeper than 40% of the continuum. 

The inferred common signature for 18 Sco is displayed in 
Eig.|^for an increasing number of spectral lines. Even when the 
full line list is considered, only a marginal detection of circu¬ 
lar polarization is possible at the radial velocity of the star. This 
is an example in which having credibility intervals for each ve¬ 
locity bin is crucial. The oscillations detected by the standard 
LSD method fall inside the uncertainty region obtained from the 
Bayesian LSD method, although in this case the mean is much 
smoother. 


5. Conclusions 

We have described a fully Bayesian least-squares deconvolution 
method which relies on a Gaussian Process prior to regularize 
the spectral shape of the common signature. Two obvious ad¬ 
vantages emerge from this Bayesian approach. Eirst, the final an¬ 
swer to the problem is in the form of a probability distribution. 
This allows the user to check whether the detected signals are 
relevant or just a consequence of fluctuations of random and/or 
systematic effects. Additionally, given that we also compute the 
covariance matrix of the ensuing probability distribution (taking 
into account possible correlations among all velocity bins), it can 
be used in inversion codes that fit the Bayesian LSD profile to 
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estimate the stellar magnetic field. This introduces a strong reg¬ 
ularization that avoids fitting spurious effects. Examples of this 
are features of the LSD profiles that are close to the noise level 
or systematic effects like interference fringes. Second, the hi¬ 
erarchical approach using a GP prior automatically adapts to the 
complexity of the spectral profile. The solution for very noisy ve¬ 
locity bins (for instance, close to the continuum) is very smooth 
and the noise is never overfitted. This automatic adaptation to the 
signal is allowed due to the inclusion of the two hyperparameters 
of the GP prior. They are easily inferred from the data. 

By taking advantage of standard linear algebra identities 
and techniques for working with sparse matrices, the Bayesian 
LSD common signature can be computed very fast in a mat¬ 
ter of a few seconds for thousands of spectral lines. The code 
to compute the Bayesian LSD profile is freely available at 
https://github.com/aasensio/pyGPLSD. 
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